Possible niche compression and individual specialization in Pacific Arctic beluga (Delphinapterus leucas) from the 19th to 20th century

Abstract Cetaceans have shown a potential to be used as sentinel species for tracking environmental change in marine ecosystems, yet our assessment of change is typically limited to recent decades and lacks ecological baselines. Using historical museum specimens, we compared community niche metrics and degree of individual dietary specialization in groups of Pacific Arctic beluga (Delphinapterus leucas) from the 1800s (n = 5) to 1900s (n = 10) using stable carbon and nitrogen isotopes drilled from teeth. Beluga occupied a broader trophic niche and demonstrated a higher degree of individual specialization in the 1800s than in the 1900s. The cause of this shift is difficult to confirm given long timescales and constraints of specimen‐based research but could indicate changes in the prey base or competition. The scale and nature of this detected shift provide perspective for continued research on these climate‐vulnerable species.

period of environmental change and human pressures (e.g, shipping, hunting, overfishing) will depend on a multitude of factors including not only the type, magnitude, and persistence of these stressors but also the biology, ecology, and behavioral responses displayed by species at the individual and population levels.
Beluga whales (Delphinapterus leucas) are a particularly valuable sentinel species for understanding the magnitude and impact of environmental change in the Arctic and sub-Arctic due to their long lives, circum-Arctic distribution, and association with sea ice (Stewart & Stewart, 1989). Additionally, they are both culturally and nutritionally important for Arctic indigenous groups. Multiple populations currently inhabit the Pacific Arctic and include the (1) Eastern Chukchi Sea, (2) Beaufort Sea, (3) Bristol Bay, (4) Cook Inlet, (5) Eastern Bering Sea, and (6) Anadyr populations (Nelson et al., 2018;O'Corry-Crowe et al., 2018). As a species, beluga have a generalist diet, with whales consuming a wide variety of prey including fish, cephalopods, and invertebrates, although Arctic cod (Boreogadis saida) is perhaps the most important prey item, especially for northern populations (Loseto et al., 2009;Quakenbush et al., 2015). The degree of individual dietary specialization in beluga is not well understood, and certain populations likely exhibit more dietary specialization than others. Variation in foraging behaviors has been documented between and even within populations depending on age and sex (Marcoux et al., 2012;Quakenbush et al., 2015).
There is currently a large degree of uncertainty as to how climate change and warming oceans are affecting beluga populations in the Pacific Arctic and how they will continue to affect beluga in the coming decades. Already, multiple pieces of evidence are suggesting that these whales are exhibiting signs of response to environmental change, including a decline in growth rates, alterations of migration routes and/or timing, and dietary switches (Bailleul et al., 2012;Choy et al., 2020;Harwood et al., 2015;Hauser et al., 2017;O'Corry-Crowe et al., 2016). Understanding how beluga have responded to changing conditions in the past is likely to improve predictions of how their populations in the Pacific Arctic might respond to current and future environmental changes. Museum specimens, including bones, teeth, hair, feathers, and baleen can be used to examine past ecosystem conditions and species behavior. Beluga teeth grow continuously over the course of an individual's life and are composed of identifiable growth layer groups (GLGs), with approximately one GLG being deposited each year (Luque et al., 2007;Nelson et al., 2018;Waugh et al., 2018). These GLGs are formed in metabolically inert dentin which remain unchanged once deposited, so a single tooth can be sequentially sampled to provide a plethora of data points over the course of a beluga's life with a temporal resolution of 1 year. Other studies have used sequential sampling of metabolically inert tissues over continuous timescales in arctic marine mammals, including beluga, narwhals, bowhead whales, ringed seals, and harp seals (Boucher et al., 2020;Dietz et al., 2021;Kershaw et al., 2021;Nelson et al., 2018;Schell, 2000).
Stable isotope analysis (SIA) has proven to be an invaluable tool for studying the biology, ecology, and physiology of living organisms, including difficult-to-study marine mammals Fry, 2006;Newsome, Martinez del Rio, et al., 2007;Smith et al., 2022;Wolf et al., 2009) as ratios of isotopes behave, at least to some degree, in predictable ways based on life history traits, environmental conditions, and species physiology. Stable isotope ratios of carbon ( 13 C/ 12 C) and nitrogen ( 15 N/ 14 N) are particularly informative for foraging ecology, as carbon is indicative of patterns of primary production and therefore foraging habitat (e.g., coastal and benthic areas are enriched in δ 13 C over offshore and pelagic areas), and nitrogen is reflective of the trophic level of consumer species and nitrogen processes at the base of the food web (Ben-David & Flaherty, 2012;Kurle & McWhorter, 2017). Patterns of δ 13 C in the ocean are dictated by primary producer community structure and growth dynamics, resource availability, and sources of dissolved inorganic carbon (DIC) as well as by terrestrial input (e.g., rivers) and, during the last century, the burning of fossil fuels by humans (known as the Suess effect; Barnes et al., 2009;Gradinger, 2009;Vega et al., 2019). δ 15 N trends are driven by the availability and dominance of certain nitrogen species in seawater, where higher levels of denitrification will lead to more enriched δ 15 N ratios than areas with higher levels of nitrogen fixation (Cline & Kaplan, 1975;Kurle & McWhorter, 2017). Additionally, and of particular importance for dietary studies, nitrogen undergoes substantial fractionation during metabolic processes, leaving tissues of a consumer more enriched in the heavier 15 N isotopes with each increase in trophic level (Brault et al., 2019;DeNiro & Epstein, 1981;Minagawa & Wada, 1984).
Isotope ratios in a consumer's tissues are reflective of diet during the period of time in which the tissue was synthesized; in this sense, δ 13 C and δ 15 N from beluga GLGs represent yearly averages of dietary information for an individual. In addition, as both δ 13 C and δ 15 N are driven by an animal's habitat and resource utilization, together these ratios can be used as a two-dimensional representation of the ndimensional niche which an individual occupies (Hutchinson, 1957;Newsome, Martinez del Rio, et al., 2007).
Many species of marine mammals exhibit high degrees of behavioral plasticity and individuality in how they occupy their ecological niche (Fleming et al., 2016;Hutchinson, 1957;Ramp et al., 2015;Samarra & Miller, 2015;Schwarz et al., 2021). A combination of pressures exerted on a consumer (e.g., intra or interspecific competition, predation, prey availability), their behavioral response to them, and other extrinsic environmental factors all act to shape the realized trophic niche at both individual and population levels (Bolnick et al., 2010;Hayden et al., 2019). For example, strong intraspecific competition among a predator species might be expected to result in a greater niche width at the population level and a higher degree of trophic variability (or dissimilarity) between individuals of the population (Bolnick et al., 2010).
A set of community niche metrics have been developed to use stable isotope ratios to quantify food web characteristics, describing the niche widths of communities as a whole as well as how groups of taxa are packed within that overall niche in carbon-nitrogen isospace (Layman, Arrington, et al., 2007). These metrics have been used for studies examining how trophic niches are affected by habitat fragmentation, novel species invasions, and man-made ecological disasters, among others (Andrades et al., 2020;Jackson et al., 2012;Layman, Quattrochi, et al., 2007). Yet these same metrics also lend themselves well to studies of high trophic level consumers, and when applied to beluga are indicative of the occupied niche width at the population level and individual dietary specialization (Layman, Arrington, et al., 2007).
The degree of dietary specialization of an individual can be calculated as the niche occupied by an individual when compared to the total niche width available to the individual (defined here as the niche occupied by the population as a whole ;Bearhop et al., 2004;Bolnick et al., 2002). An individual is considered to have a specialist diet when its realized niche width is relatively small compared to its total available niche width, while an individual occupying a larger niche width relative to its total available niche width would be considered to have a more generalist diet (Araújo et al., 2007;Newsome et al., 2009). This conceptualization of individual dietary specialization has been applied to a number of taxa, including southern sea otters (Enhydra lutris nereis), southern elephant seals (Mirounga leonina), bull sharks (Carcharhinus leucas), tiger sharks (Galeocerdo cuvier), and leptodactylid frogs, by using stable isotope ratios of δ 13 C and δ 15 N taken from tissues (Araújo et al., 2007;Hückstädt et al., 2012;Matich et al., 2011;Newsome et al., 2009). The degree of individual dietary specialization within a species can have important ramifications regarding how that species might persist in the face of environmental change and other anthropogenic stressors.
In this paper, we aim to describe beluga foraging ecology in the Pacific Arctic from the 1800s to the 1900s using δ 13 C and δ 15 N values from individual tooth GLGs in an attempt to understand how populations may have changed in the past and consequently how they may continue to adapt in an uncertain future. Specifically, our goals are to (1) compare community isotopic niche metrics between beluga of both centuries and (2) determine the level of individual specialization in beluga by comparing the within-individual isotope variability from samples taken from each GLG to the within-century variation across all individuals. This study will help place modern studies of beluga into context by providing historical baseline isotope data and may allow for better-informed management of the system and species in the future.

| Sample collection
To ensure a broad temporal coverage of beluga tooth isotope data from the Pacific Arctic, teeth were sourced from multiple locations throughout the Pacific Arctic and from two distinct time periods.
Teeth with GLGs ranging in date from the latter half of the 19th century (approximately 1845-1898; n = 5) and the mid-20th century (1940-1983; n = 10) were loaned with permission for destructive sampling from the Smithsonian National Museum of Natural History in Washington, DC and the University of Alaska's Museum of the North (UAM) in Fairbanks, Alaska. Specific localities of harvested, captured, or found beluga, as well as temporal coverage of data from these individuals, can be found in Figure 1 and Table 1. Specifics regarding original collection methods are unknown, although it is very likely that 1800s samples came from early scientific expeditions and 1900s samples were sourced from indigenous harvests of whales during migrations. While localities were included for all the specimens, further details on population assignment were unfortunately not available. Sex was unknown for all but one beluga from the 1800s; a young female which was 8 years old at death. Out of the 10 whales from the 1900s, there were four males, four females, and two individuals of unknown sex.

| Stable isotope analysis
Upon receipt from museum collections, teeth were cut in half longitudinally to reveal GLGs. Powder was drilled from every GLG when possible, providing a temporal resolution of 1 year. For some teeth with extremely narrow GLGs two, three, or in one single case four layers were combined, resulting in a total of 304 drilled samples. When multiple GLGs were combined, a "mean year" value was given for that sample by taking the mean of included years and rounding up to the nearest year. Drilling was done with either a computer-guided micromilling system (Carpenter Microsystems

| Statistical analyses
All statistical analyses were conducted using R (version 4.0.2; R Core Team, 2020). Wilcoxon tests were used to compare δ 13 C and δ 15 N between males and females from the 1900s beluga, with mean values per individual (n = 4 females and n = 4 males). Of the 1800s beluga, only one had known sex, a female, so no analyses were made for individuals in this group. Due to the opportunistic nature of acquiring museum specimens for research, individual beluga were harvested over a large geographic area with differential spatial and temporal coverage across the two time periods (1800s and 1900s; see Figure 1).
Given this uneven sampling across centuries and the nature of isotope signatures, particularly carbon, at the base of the food web to vary predictably through space and time (Goericke & Fry, 1994;Kurle & McWhorter, 2017), we first wanted to explore the potential scope and relative influence of both location and year on the beluga isotope signatures. A pairwise comparison of all sampled points was made, where isotope ratios of every single GLG were compared with every other GLG in the dataset. This resulted in a new dataset with the change of δ 13 C and δ 15 N (∆ δ 13 C and ∆ δ 15 N) between each pair, as well as the geographic distance between paired samples (zero for GLGs from the same whale or multiple whales taken from the same location) calculated using the package swfscMisc (Archer, 2022) and the difference in time (years) between paired samples. The effect that geographic distance and difference in time had on the change in δ 13 C or δ 15 N was addressed using linear regressions. F I G U R E 1 (a) Map of the study area and collection locations of teeth. Gray circles represent teeth from the 1800s and green circles represent teeth from the 1900s. (b) Temporal coverage of each individual beluga.
Community niche metrics, hereafter referred to as Layman metrics, were calculated using the package SIBER (Jackson et al., 2011;Layman, Arrington, et al., 2007). All six of the following metrics were calculated per model iteration, and the modes and 95% credible intervals, the latter of which were calculated as the highest density region, of the resulting posteriors are reported here. Total area (TA) represents the minimum area in isospace (the two-dimensional space with δ 13 C as the x-axis and δ 15 N as the y-axis) which a polygon occupies while including centroid locations for each individual of a community (in this case century). Range in δ 13 C (CR) and δ 15 N (NR) is the one-dimensional difference between maximum and minimum values for each isotopic ratio and is similar to TA in that it is reflective of trophic niche width in a community. The mean distance to the centroid (CD) is the mean distance of each individual beluga to the mean of all beluga in the given community (century). Last, the mean distance to the nearest neighbor (MNND) and standard deviation distance to the nearest neighbor (SDNND) are the mean and standard deviation, respectively, of each beluga to its nearest neighbor in isospace. Bayesian models with vague priors were fitted to each individual beluga to calculate these metrics as well as stan- To calculate the degree of individual specialization in beluga the following formula was used: S = WIC/(BIC + WIC), where S is the specialization index, WIC is the within-individual component, BIC is the between-individual component, and the sum of BIC and WIC is the total niche width of the population (Bearhop et al., 2004;Hückstädt et al., 2012;Newsome et al., 2009). The smaller the value of S, the more specialized an individual is considered in relation to the total niche width available to that individual. Previous work has considered an S ≤ 0.2 to be an extreme specialist whereas an S ≥ 0.5 would indicate a generalist individual (Hückstädt et al., 2012). The WIC and BIC are calculated as the standard deviation of a set of points-here the WIC is the standard deviation of all sampled GLGs in an individual whale and the BIC is the standard deviation of all sampled GLGs of all beluga in the respective time period (1800s vs 1900s). It is important to note that the degree of individual specialization can be affected not only by changing an individual's WIC, but also by an expansion or contraction of the population's total niche width.

| RE SULTS
Mean δ 13 C for all beluga was −13.72‰ and mean δ 15 N was 19.18‰ ( Figure 2). When broken down by century, 1800s beluga had a mean TA B L E 1 Metadata for all 15 belugas included in this study. Catalog numbers are as listed in the respective museum (USNM, Smithsonian Institution's National Museum of Natural History & UAM, University of Alaska Museum of the North). Samples drilled do not necessarily match the number of years of coverage, as some GLGs were combined together into one sample.  (Woodring, 1958). beluga predicted isotope ratios from 1800s beluga revealed some distinction between centuries for both ratios; for δ 13 C the mode of 1800s-1900s δ 13 C and δ 15 N was 0.129‰ (95% credible interval −0.184‰ to 0.441‰) and −0.350‰ (−0.669‰ to −0.039‰) respectively. The probability that δ 13 C and δ 15 N from 1800s beluga were greater than that of 1900 beluga, calculated as the percent of the posterior ≥0 was 79.5% for δ 13 C and 1.4% for δ 15 N. Conversely, the probability that 1800s δ 13 C and δ 15 N was lower than those from the 1900s was 20.5% and 98.6%, respectively.

| Pairwise comparison analysis
With a total of 304 sampled GLGs, the pairwise comparison resulted in a set of 46,056 paired observations. Linear models for both δ 13 C and δ 15 N showed a statistically significant relationship between geographic distance and time between samples and resulting isotope ratios (p < .001). However, the models estimated that on average, 1 km is related to a change in δ 13 C of 0. for δ 15 N. Thus, differences in δ 13 C are greater over distance and δ 15 N changes more over time, as expected, yet both vary by levels comparable to differences among individuals sampled in the same year and location, suggesting neither variable on its own is predictive in a linear manner nor biologically meaningful in the differences observed across the beluga GLG results.

| Layman metrics and standard ellipse areas
Gelman-Rubin diagnostics for each model indicated strong convergence of MCMC chains and confidence in the resulting posterior draws. Additionally, autocorrelation along chains was found to be low (≤0.01) when checked at lags of 100, 500, 1000, and 5000. All six Layman metrics were greater for 1800s beluga than 1900s beluga. Modes and Bayesian credible intervals are reported in Table 2.
The greatest discrepancies between the two centuries were found in TA (the total area in carbon-nitrogen isospace) and the range in δ 15 N, the latter of which contributed more to TA than did the range

| Individual specialization
The degree of individual specialization in Pacific Arctic beluga, defined as each individual's observed occupied niche compared to the population's total available niche, was on average greater in the 1800s than in the 1900s (Table 3, Figure 4)

| DISCUSS ION
The samples analyzed here suggest that at the community level, Teeth included in our study from Kotzebue Sound include data from 1983 at the very latest, conveniently timed prior to the mysterious population decline, and so likely belonged to beluga from the Kotzebue Sound stock only. Thus, despite collection locations for 1800s beluga covering a wider geographic area than the 1900s, we still might expect a similar level of stock diversity included in each F I G U R E 4 Individual specialization index calculated from carbon and nitrogen ratios. The lower the value, the more specialized the individual whale. Overall, individual beluga appear to have utilized more of their total available resource pool during the 1900s as the S index, and level of generalization of diet, increased.
temporal group. It is also worth noting the individuals with collection locations of Point Hope, AK, and northwards most likely belong to the Chukchi and Beaufort Sea stocks. These beluga undergo some of the longest migrations of all Pacific Arctic stocks, and as a result, might be expected to have larger standard ellipse areas due to the inclusion of prey items from drastically different summer and winter areas. Yet the individual with the most northern collection location of Nuvuk (Point Barrow; catalog ID 305071 from the 1800s group) had the smallest SEA of all included whales, regardless of century ( Figure 3). Of the remaining 1900s beluga, there did not appear to be any patterns of SEA between individuals collected from Point Hope and Point Lay versus those from Kotzebue Sound (Figure 3). Thus, despite likely having considerably larger home ranges, the northernmost beluga examined in our limited set of specimens did not necessarily show a more generalist diet than the rest as indicated by stable isotope ratios.
We found female beluga to be significantly depleted in δ 15 N when compared to males, though we would not expect this to create major biases within temporal groups. For the 1900s beluga, there were approximately the same number of males and females. For the 1800s beluga, only one whale was known to be female while the rest were unknown, preventing any inferences from being made. However, this one female from the 1800s group had among the highest δ 15 N of all included whales (Figure 2; catalog ID 22207). Finally, age at death between belugas from the two time periods was not even; three individuals from the 1800s group were considerably younger at death than all beluga from the 1900s group. In female southern elephant seals, dietary specialization has been shown to increase with age (although not necessarily linearly) as individuals home in on foraging strategies (Hückstädt et al., 2012). Our analyses do not allow us to examine the degree of dietary specialization with age, although it should be noted that there were no discernable patterns in individual standard ellipse areas by age at death (Figure 3). Nonetheless, it is possible that uneven distributions of males and females, as well as the age at death of individuals, between both time periods could affect our results.
Our findings of a compressed niche width and higher levels of dietary generalization in 1900s beluga could also be explained by a multitude of scenarios concerning competition for resources, changes in beluga prey diversity and availability (ecological opportunity), large-scale shifts in biological communities (e.g., invasion or northward expansion of novel species), or predation. Given the limits of relying exclusively on bulk stable isotope ratios of carbon and nitrogen taken from teeth of historic samples and uncertainties regarding the Pacific Arctic marine environment from as far back as the mid-1800s, it is impossible to rule out any one of these scenarios in favor of another. As ratios of δ 13 C and δ 15 N are assimilated into the tissues of a consumer from its total diet during the period of tissue synthesis, any change in niche width of a consumer is ultimately due to a change in the diversity of prey consumed (Layman, Quattrochi, et al., 2007). This could be the result of consumers foraging on a larger or smaller subset of prey or an expansion or contraction in the trophic diversity of the food web as a whole; Layman, Quattrochi, et al. (2007) found that a less diverse set of basal resources was a principal driver in the collapse of the niche width of a top predator in fragmented aquatic ecosystems. While historical isotope data from the base of the food web are not available, future research could incorporate compoundspecific stable isotope analysis to allow further inferences regarding the occurrence, magnitude, or effect that a change in the baseline δ 13 C & δ 15 N of the Pacific Arctic marine environment may have had on the isotope ratios of our beluga teeth.
Factors such as prey availability or diversity, competition between and within species, and predation risk (increasingly an issue due to northward expansion of orcas) can independently and synergistically affect consumer trophic dynamics and would also presumably affect the degree of individual specialization among individuals (reviewed in Araújo et al., 2011). For example, competition for resources between and within species can either increase or decrease the trophic niche width of a consumer. Optimal foraging theory postulates that when intraspecific competition for a limited amount of resources is high, certain individuals in a population will be pressured to begin selecting for suboptimal prey types, thereby increasing the overall niche width at the population level and increasing dietary specialization and differentiation at the individual level (Kernaléguen et al., 2015;Ratcliffe et al., 2018;Tinker et al., 2012;Yurkowski et al., 2016). If the effects of intraspecific competition and prey availability as they relate to optimal foraging theory are at play here, beluga may have shown a compressed niche and higher degree of dietary overlap due to reduced competition for resources over time. In this hypothetical scenario, a release from high levels of intraspecific competition (either from less beluga or more prey) would facilitate more individuals of a population being able to select for higher quality prey and a subsequent decrease in populationlevel niche width as well as a decrease in the degree of individual specialization.
Any changes in the Pacific Arctic marine food web of a sufficient magnitude to affect the isotopic profiles of beluga prey, or simply the addition or removal of prey species, will affect the isotopic ratios found in these beluga teeth. Warming and increasingly  (Yurkowski et al., 2018). The authors argued that this decrease in trophic diversity was primarily driven by a reduction in the diversity of basal resources (similar to Layman, Quattrochi, et al., 2007) due to predators eating novel prey (e.g., capelin) and switching largely to a pelagic energy pathway.
Capelin have been documented as far north as the Beaufort Sea as early as the 1960s; this encroachment is quite recent relative to the timeline covered by our 1900s samples and may not be as relevant to our results (Choy et al., 2020;McNicholl et al., 2017). Still, all but one beluga from our 1900s community had a date of death of approximately 1977-1983 so it is possible that the encroachment of capelin, and/or similarly novel species moving north as oceans warm, could have been at play here. Finally, the Pacific Arctic experiences strong, decadal-scale regime shifts due largely to climatic drivers (e.g., Pacific Decadal Oscillation, North Pacific Gyre Oscillation, and El Niño-Southern Oscillation, among others), where species assemblages are significantly altered (Brodeur et al., 2008;Conners et al., 2002;Ebbesmeyer et al., 1991;Hare & Mantua, 2000;Litzow et al., 2014;Niebauer, 1998). These regime shifts likely have a profound effect on the biology of marine mammal predators which may be reflected in isotope ratios (see York et al., 2008). Additionally, numerous studies have detected long-term changes in primary production and carrying capacity in the region using time series data from marine mammal tissues (Newsome, Etnier, et al., 2007;Schell, 2000).
While recent studies show significant evidence for contemporary shifts in food webs and predator foraging behavior (Choy et al., 2020;Marcoux et al., 2012;Watt et al., 2016), the data presented here suggest that diet shifts may also have occurred prior to the most significant climate-driven changes in the region. It seems possible that over the 138-year timeframe of the present study, there may have been numerous shifts in the prey landscape for beluga that would have resulted in isotopic shifts or changes in niche width. The scale and nature of this trophic shift may provide helpful context for contemporary observations that aim to understand the increasing impact of climate change on beluga whale diet, health, and behavior. On one hand, the shifts documented here suggest beluga have already been quite plastic in their trophic behavior, both in terms of individual specialization and community niche breadth. While that suggests they may continue to adapt in the face of current climate-driven changes, the compression of niche width may also indicate an increasingly limited foraging landscape for the species which may decrease resiliency. As one of just a few endemic Arctic cetaceans and a culturally and nutritionally significant resource for Alaska Natives, further work on the foraging behavior of beluga and their Pacific Arctic food web should be prioritized.

Smithsonian Institution's James Smithson Fellowship and National
Science Foundation grant #1927742.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.3bk3j 9kqp.